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Abstract 

We further develop a thermal LB model for multiphase flows. In the improved model, we 
propose to use the windowed FFT and its inverse to calculate both the convection term and 
external force term. By using the new scheme, Gibbs oscillations can be damped effectively in 
unsmooth regions while the high resolution feature of the spectral method can be retained in 
smooth regions. As a result, spatiotemporal discretization errors are decreased dramatically 
and the conservation of total energy is much better preserved. A direct consequence of 
the improvements is that the unphysical spurious velocities at the interfacial regions can be 
damped to neglectable scale. With the new model, the phase diagram of the liquid-vapor 
system obtained from simulation is more consistent with that from theoretical calculation. 
Very sharp interfaces can be achieved. The accuracy of simulation results is also verified by 
the Laplace law. The high resolution, together with the low complexity of the FFT, endows 
the proposed method with considerable potential, for studying a wide class of problems in 
the field of multiphase flows and for solving other partial differential equations. 
Keywords: Lattice Boltzmann method; spurious velocities; liquid-vapor systems; 
windowed FFT 
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1. Introduction 



During the past two decades, the lattice Boltzmann C 



(B) method has been developed 



rapidly and has been successfully applied to various fields [1] , ranging from magnetohydrody- 
namics {2, 3, 4], to flows of suspensions jsl, 6], flows throu gh p orous media (3, IJ compressible 

16] . etc. Apart from the fields 



fluid dynamics [9 



10 



11 



12 



wave propagation 



15 



isted above, this versatile method is particularly promising in the area of multiphase systems 



17 



18 



19 



20 



21 



22 



23|. This is partly owing to its intrinsic kinetic nature, which makes 
the inter-particle interactions (IPI) be incorporated easily and flexibly, and, in fact, the IPI 
is the underlying microscopic physical reason for phase separation and interfacial tension in 
multiphase systems. 

So far, several LB multiphase models have been proposed. Among them, the four well- 



known models are the chromodynamic model by Gunstensen et al. [17|], the pseudo-potential 



model by Shan and Chen (SC) 
model by He, Chen, and Zhang 



the free-energy model by Swift et al. [19|, and the HCZ 
201 ] . The chromodynamic model is developed from the two- 



component Lattice Gas Automata (LGA) model originally proposed by Rothman and Keller 



24J. In this model, the red and blue colored particles are employed to represent two different 



fluids. Phase separation is achieved through controlling the IPI based on the color gradient. 
Similar to the treatment in molecular dynamics (MD), in SC model, non-local interactions 
between particles at neighboring lattice sites are incorporated. The interactions determine 
the form of the equation of state (EOS). Phase separation or mixing is governed by the 
mechanical instability when the sign of the IPI is properly chosen. In the free-energy model, 
besides the mass and momentum conservation constraints, additional ones are imposed on 
the equilibrium distribution function, which makes the pressure tensor consistent with that 
of the free-energy functional of inhomogeneous fluids. In the HCZ model, two distribution 
functions are used. The first one is used to compute the pressure and the velocity fields. The 
other one is used to track interfaces between different phases. Molecular interactions, such 
as the molecular exclusion volume effect and the intermolecular attraction, are incorporated 
to simulate phase separation and interfacial dynamics. 

The aforementioned models have been successfully applied to a wide variety of multiphase 
and/or multicomponent flow problems, including drop breakup 25, 26], drop collisions 2?|, 



wetting 



2n 



29], contact line motion 



tion and phase ordering phenomena 



30 



19 



31 



21 



, chemically reactive fluids 



23 



321 ]. phase separa- 



331 ] . hydrodynamic instability (3J,]35|, etc 



However, despite this, the current LB versions for multiphase flows are still subjected to, 
at least, one of the following constrains (i) the isothermal constraint (i.e., the deficiency of 
temperature dynamic), (ii) the limited density ratio and temperature range, (iii) the spurious 
velocities. This paper addresses mainly the last restriction and the total energy conservation 
in practical simulations. 

Spurious velocities extensively exist in simulations of the liquid-vapor system and reach 
their maxima at the interfacial regions, indicating deviation from the real physics of a fluid 
system. Reducing and eliminating the unphysical velocities are of great importance to the 
simulations of multiphase flows. Firstly, large spurious velocities will lead to numerical 
instability. Secondly, the local velocities are small during phase separation and coarsening. 
If the spurious currents are too large, then we may not be able to separate the spurious 
currents from the real local flows, which is especially true in the case of phase separation 
with high viscosity. Thirdly, for a thermal multiphase system, accurate flow velocities are 



36]. 



required in order to obtain an accurate temperature field 

In dealing with this issue, extensive efforts have been made during the past years. Wagner 
[37I pointed out that the origin of the spurious currents is due to the incompatibility be- 
tween the discretizations of driving forces for the order parameter and momentum equation. 
Therefore, he suggested to cure the spurious velocities by removing the nonideal terms from 



23. 



38 



the pressure tensor and introducing them as a body force. Sofonea and Cristea et al. 
presented a finite difference LB (FDLB) approach and proposed two ways to eliminate the 
unwanted currents. In the first way, a high-accuracy numerical scheme, the flux limiter 
method is employed to calculate the convection term of the LB equation. In the second way, 
a correction force term is introduced to the LB equation that cancels the spurious velocities 
and allows to recover the mass equation correctly. Shan 39j and Succi et al. 22J showed 
that the origin of the spurious currents are due to the insufficient isotropy in the calculation 
of density gradient. Therefore, using the information of the density field on an extended 
neighboring x + e» of a given site to construct high order isotropic difference operators, is 



the key for the correct discretization of spatial derivatives and taming the spurious currents 



in the interface. Yuan et al. [36] demonstrated that smaller parasitical velocities and higher 
density ratio can be achieved using more realistic EOS in a single-component multiphase 
LB model. Lee and Fischer 40] reported that the use of the potential form of the surface 
tension and the isotropic FD scheme can eliminate parasitic currents to round-off. Seta and 



Okui 



41] composed a more accurate fourth order scheme to calculate the derivatives in the 



pressure tensor. This convenient approach reduces the amplitude of spurious velocities to 
about one half of that from the second order scheme. Pooley and Furtado 42j analyzed 
the causes of spurious velocities in a free-energy LB model and provided two improvements. 
First, by making a suitable choice of the equilibrium distribution and using the nine-point 
stencils (NPS) scheme to calculate derivatives, the magnitude of spurious velocities can be 
decreased by an order. Moreover, a momentum conserving force is presented to further 
reduce the spurious velocities. Yeomans et al. [43] identified two sources of the spurious 
velocities, the long range effects and the bounce-back boundary conditions, when a single 
relaxation time (SRT) LB algorithm is used to solve the hydrodynamic equations of a binary 
fluids. Aiming to reduce the unwanted velocities, they proposed a revised LB method based 
on a multiple-relaxation-time (MRT) algorithm. 

In this work, we present a thermal LB model for simulating thermal liquid- vapor system 
with neglectable spurious velocities. This model is a further development of the one originally 
proposed by Watari and Tsutahara (WT) [tj and then developed by Gonnella, Lamura and 



Sofonea (GLS) [44] . The original WT model works only for ideal gas. GLS introduced an 
appropriate IPI force term to describe van der Waals (VDW) fluids. Here we introduce a 
windowed FFT (WFFT) scheme to calculate the convection term and the force term. The 
improved model is convenient to compromise the high accuracy and stability. With the 
new model, non-conservation problem of total energy due to spatiotemporal discretizations 
is much better controlled and spurious currents in equilibrium interfaces are significantly 
damped. 

The rest of the paper is structured as follows. In the next section the thermal LB 
models for ideal gas and for VDW fluids are briefly reviewed. In section III we illustrate 
the necessity of the further development and detail the usage of the WFFT scheme and its 
inverse. Comparisons and analysis of numerical results from different schemes are presented 



in section IV, where we will show how the spurious velocities around linear and curved 
interfaces can be reduced by the new model. Finally, in section V, we summarize the results 
and suggest directions for future research. 



2. The model 

2.1. Original WT model for ideal gas system 

The thermal multiphase model is developed from the thermal LB model, originally, pro- 
posed by WT, which is based on a multispeed approach. In this approach, additional speeds 
are required and higher order velocity terms are included in the equilibrium distribution 
function to obtain the macroscopic temperature field. 

WT model uses the following discrete-velocity-model (DVM) which involves a set of 33 
nondimensionalized velocities 

i — i i — 1 

v = 0, v fci = Ufc [cos(— ^— 7r),sin(— ^— tt)], (1) 

where subscript k = 1,...,4 indicates the fc-th group of the particle velocities whose speed 
is Vk and i = 1,...,8 indicates the direction of particle's speed. In our simulations we set 
vi = 1.00, v 2 = 1.90, v 3 = 2.90, and v± = 4.30. 

The distribution function f^, discrete in space and time, evolves according to a SRT 
Boltzmann equation 

9fki dfki _ 1 r f req-i / n \ 

~df Wki ' ~~dr~ ~ ~t ^ ~ ' ^ ' 

where f%f, r, and r denote the local equilibrium distribution function, the spatial coordinate, 
and the relaxation time, respectively, f%jf is expressed as a series expansion in the local 
velocity 

req P ( V li\ ( U 2 - 2UV ki 
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where the weight factors are 

Fk = « - ™ - <m - «w [48T 2 4 - «^ + ^ + <3)^ 

+«i< 2 +< 2 < 3 +< 3 <i)^ 2 - <A^ T], (4) 



F = l-8(Fi+F 2 +F3+F 4 ), (5) 



with 



fc + Z, /c + Z<4; 
{£; + /}=<( " Ze {1,2,3}. (6) 

fc + Z-4, k + l>A; 

Hydrodynamic quantities, such as density, velocity, and temperature are determined from 
the following moments 

P = £/*1, (7) 

ki 

P U = ^Z V kifki, (8) 

ki 

pr = X^(v w -u) a /£. (9) 

The combination of the above DVM and the general FD scheme with first-order forward 
in time and second-order upwinding in space composes the original FDLB model by WT. 
In the FDLB model, particle velocities are independent from the lattice configuration. As a 
result, higher-order numerical schemes can be used to reduce the numerical viscosity and to 
enhance the stability of the model. This is of great importance to LB simulations, especially 
in phase separation studies, where long lasting simulations are needed to establish the growth 
properties. 

2.2. GLS model for multiphase system 

WT model can be applied to compressible flows with small Mach number and the revised 
version [2] extends it to compressible flows with high Mach number due to better numerical 
stability. Nevertheless, neither the original one nor the improved one has the ability to 
describe multiphase flows, since both models lead to the ideal EOS only, which do not 
support thermodynamical two-phase state. 



Fortunately, by incorporating a forcing term, the improved model can be applied to 
thermal liquid-vapor systems. Compared to isothermal models, the variable temperature 
that the GLS model can be implemented is of great importance, since thermal effects are 
ubiquitous and sometimes dominant in an important class of flows 45[. Examples are referred 
to boiling 46[, distillation, as well as the dynamics of phase separation 47|, |48j, where the 
freedom in temperature limits the rate of phase separation and induces different rheological 
and morphological behaviors. Dynamic effects of temperature can not be considered in 
isothermal models, therefore, most studies have been restricted to either isothermal systems 
or the systems where effects of temperature dynamics are negligible. 

The forcing term introduced by GLS is added into the right-hand side (RHS) of Eq. (2) 



dfki dfki 
ot or 

where takes the following form, 



[fki fki) Ikii 

T 



(10) 



■[A + B a (v kia - u a ) + (C + C q )(v kia - u a ) 2 ]f^. 



Iki in Eq. (10) is introduced to control the equilibrium properties of the liquid-vapor systems 
and allows to recover the following equations for VDW fluids 44|, 



where 



d t p + d a (pu a ) = 0, 

d t {pu a ) + dp(pu a ui3 + n a/3 - a a p) = 0, 
d t e T + d a [e T u a + (U a/3 - a a p)u p - K T d a T] = 0, 

na/3 = P W 5 a/ 3 + A a p, 



(12) 

(13) 
(14) 

(15) 
(16) 



and 



e T = P T-9p 2 /8 + K\Vp\ 2 /2 + pu 2 /2, (17) 

represent the non-viscous stress, the dissipative tensor, and the total energy density, respec- 
tively, kt, rj, and ( are heat conductivity, shear, and bulk viscosities. P w and A Q/ 3 in Eq. 



([15]) are the VDW EOS and the contribution of density gradient to pressure tensor, which 
have the following expressions 

P w = 3pT/(3 - p) - 9p 2 /8, (18) 

A Q/3 = Md aP d pP - [pTd,pd,(M/T)]5 a p - M(pV 2 p + | Vp| 2 /2)6 af) . (19) 

The expression M = K + HT allows a dependence of the surface tension on temperature, 
where K is the surface tension coefficient and if is a constant. 

In order to recover Eqs. f)T2]) - ( TT^|) . five constraints are imposed on the forcing term, which 
make coefficients in Eq. ( TTTj) as the following form 

A = -2{C + C q )T, (20) 

B a = -^[d a (P w - pT) + dpKe - d a ((d lUl )\, (21) 

C = Yjv{( pW ~ P T )d~i u i + ^apd a u p - (C<9 7 M 7 )9 Q w a 

9 1 
+-p 2 a 7 w 7 + K[--(d 1 p)(d 1 p)(d a u a ) 

O Z 

-p(d^p)(d y d a u a ) - (d y p)(d y u a )(d a p)]}, (22) 

C q = -^d a [2qpT(d a T)}. (23) 

It is worth noting that in this model the Prandtl number Pr = T)/kt — t/2(t — q) can be 
changed by adjusting the parameter q in the term C q . 

3. Thermal LB model based on the WFFT approach 

In this section, we present our contribution to the thermal multiphase LB model: spatial 
derivatives in the convection term \ki • dfu/dv and in the forcing term J^j, are calculated 
via the WFFT algorithm and its inverse. 

To illustrate the necessity, we present simulation results for a thermal phase separation 
process by various numerical schemes. Here the time derivative is calculated using the 
first-order forward Euler FD scheme. The spatial derivatives in 1^ are calculated using the 
second-order central difference scheme. Spatial derivatives in the convection term Vki-dfki/dr 
are calculated using various schemes listed as follows: 



3.1. second-order central difference scheme 

Let J — 1, J and J + 1 be three successive nodes of the one dimensional lattice. Using 
the second-order central difference scheme to discretize the convention term, Eq. ([TO]) can 
be rewritten in a conservative form 

en+1 _ en C kia / fn fn \ ~_f fn eeq,n\ , T n a, (1A\ 

Jki,J — Jki,J ^~\Jki,J+l ~ Jki,J-l) ~{Jki,J ~ Jki,j) 1 ki,J lAl , 

where At and c^a = v^At/ Ar a are the time step and the Courant-Friedrichs-Levy (CFL) 
number. 

3.2. Lax-Wendroff (LW) scheme 

Compared with the second-order central difference scheme, the LW scheme contributes 
a dissipation term, which is in favor of the numerical stability. Then, by using this scheme, 
Eq. ([TO]) can be formulated as 

fn+1 m _ Ckia i p n p n \ ^kia I en _ ry en , en \ 

J ki,J ~ Jki,J 2 Jki,J-l) ~r \J ki,J+l Z J ki,J "r Jki,J-X) 

-^(fki,j-f e k!;j) + I ki,J At - ( 25 ) 

3.3. Modified-LW (MLW) scheme 

As we know, the LW scheme is very dissipative and has a strong " smoothing effect". 
Obviously, it is not favorable to recover the sharp interface in the multiphase system. To 
further improve the numerical accuracy, the modified partial differential equation (MPDE) 
remainder after discretizing with Eq. ( )25l) is derived 49 ] 



p _ ~ C L) a -2 9 3 f VkiaCkioQ- ~ 4ia) A -3 & j / 9R n 

It is clear that the first and the second terms in the RHS of Eq. (j26l) correspond to the third- 
order dispersion error R% and the fourth-order dissipation error R^, respectively. Therefore, 
we can add the dispersion term into the RHS of Eq. ([TO]) to compensate the dispersion error 

dfki 9 ^ f f _ f eq \ J- T 

1 ^ki-~Q^Jki — ~\Jki J ki ) ' * ki 

Vkiaj^- - C 2 kia ) 2 3 f 
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J-l 



Figure 1: (Color online) Characteristic line on the square lattice for the direction i = 2. 

Furthermore, we can add the dissipation term into RHS of Eq. ( 1271) . Using the 2nd-CD 
scheme to discrete -R3 and -R4 gives 



C-kia ( 1 



'kia) 



and 



Ra 



'kia 



1 



12 



■'kia) 



( -f n — 9 -f n 

Uki,J+2 ~ Z Jki 



i,J+l 



VI, 



fki,J-2)i 



(fk 



(28) 



(29) 



lki,J+2 ^fki,J+l + 6fki,J ^fki,J-l + fki,J-2J- 

The bars above R% and R4 indicate that they are discretized. If only _R 3 is added into the 
RHS of Eq. (125]) . for convenience of description, we refer to this scheme as MLW1. If both _R 3 
and R4 are added into the RHS of Eq. (l2"5j) . then a more accurate LB equation is obtained, 
and we refer to this scheme as MLW2. 



3.4- Flux limiter (FL) scheme 

The FL scheme has been widely employed by Sofonea et al. 



23 



to reduce the spurious 



velocities and to improve the numerical stability in liquid- vapor systems. Figure 1 shows the 
characteristic line on the square LB lattice for direction i = 2. When using this approach to 
compute the convective term along the characteristic line, Eq. fflOl) becomes 



f n+X 
JkiJ 



rn v k^t . ^ n ^ 

Jki,J A-Ar ki ' J+1 / 2 ki,J-l/21 

-;(/^-0 At + 4V A *> 



with 



At 



1, i €{1,3,5,7}; 
y/2, ie {2,4,6,8}. 



(30) 
(31) 

(32) 



F n and F n 

r ki,J+l/2 dI1U r ki,J-l/2 



in Eq. fl30|) are two fluxes, which are defined as 



ki,J+l/2 



fki,J + 2 (1 



VkAt 

AAr a 



)[fki,J+l fki,Jl' t / } {^ki,j) 



(33) 
(34) 



r ki,J-l/2 ~ r ki,{J-l)+l/2i 

where the flux limiter ip{6 kiJ ) is expressed as a smooth function 

/n « 

/jn _ j fei, J J ki,J—l 
U ki,J ~ ~FH _ fn ' 

Jki,J+l J ki,J 

In particular, if 1^(0^ j) = 0, it corresponds to the first-order upwind scheme and ip{0 r k 
to the LW scheme. A wide choice of flux limiters can work well with LB models. In this 
work, we will use the monitorized central difference (MCD) FL, which is most widely used 
by Sofonea et al. 



'kiJ 



(35) 
= 1 



0. 



Otj < o, 



< Btj < 1/3, 



ki,J) 



(36) 



(l + ^,j)/2, l/3<^<3, 

2, 3 < e^j. 

3.5. NPS scheme 

Recently, a new scheme, named the NPS scheme, has been widely used to calculate the 
spatial derivatives by many scholars so as to ensure higher isotropy and to reduce spurious 



velocities 
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5l|. The general choice of stencils for calculating the derivatives 



and the laplacian are 
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Figure 2: (Color online) Variations of total energy Aex(t) = er(t) — er(0) for a phase separating process 
with various numerical schemes. 

with 2 A + 4B = 1 and E + 2F = 1 to keep consistency between the continuous and discrete 
operators. The bars above d x and V 2 represent that they are discrete operators. The central 
entry denotes the lattice node at which the derivative is calculated, and the remaining entries 
are the eight neighbor nodes around the central one. B and F are two free parameters that 
are chosen to minimize the spurious velocities. A large amount of numerical tests indicate 
that the best choice is B — 1/12 and F — 1/6 in the GLS model. In our simulations, both 
the convection term and the forcing term are calculated by this way. 

Next, we conduct simulations of a thermal phase separation process with numerical 
schemes listed above. Initial conditions of our test are chosen as 



where A is a random density with an amplitude 0.001 and can be regarded as incipient 
nuclei in the density field. Periodical boundary conditions (PBC) are imposed on a square 
lattice with N x = N y = 128. Unless otherwise stated, the remaining parameters are Ax = 
Ay = 1/256, At = 10~ 5 , r = 10" 2 , K = 10~ 5 , H = 0, ( = 0, q = -0.004, throughout 
our simulations. Figure 2 shows the variations of total energy Ae^(t) = ex(t) — er(0) for 
the phase separating process with various numerical schemes. The legend in each case is 
composed of two parts, 'A'+'B', where 'A' is 'CD', 'LW\ 'MLW1', 'MLW2', 'FL' and 'NPS' 
and it shows the scheme to calculate the convection term; 'B' is 'CD' and 'NPS' and it 



(p, T, u, «) = (! + A,0.85, 0.0, 0.0), 



(39) 



shows the scheme to calculate the forcing term. Figure 2 demonstrates that the total energy 
density ey(t) is not conservative in simulations even though it is in theoretical analysis. 
Further survey of these results indicates that the derivation Aerit) decreases by increasing 
the accuracy of scheme. Therefore, we conclude that the non-conservation of total energy is 
mainly due to the spatial discretization errors. 

3.6. The WFFT approach 

To overcome the problem of energy non-conservation, a new algorithm based on WFFT 
is proposed. This approach is especially powerful for periodic system and also provides 
spatial spectral information on hydrodynamic quantities. Moreover, with this approach, 
higher-order derivatives and fractional-order derivatives can be computed in a convenient 
way. 

For the sake of clarity, we start with the definition of Fourier transform of f(xj) 

N-l 

f(k) = AxJ2f(^)e~ ikx \ (40) 

3=0 

and its inverse 

N/2-1 
n=-N/2 

where k is the module of wave vector k, i is the imaginary unit, and f(k) stands for the 
Fourier transform of a spatial function f{x). In Eq. (j4Tj) . k=2irn/L, and L = NAx is the 
length of the system divided into N equal segments. The above two equations are exactly 
correct when N is infinitely large or Ax is infinitely small. A general theorem of derivative 



based on FFT states that 
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54) 



/'(*) = ik x /(*), (42) 



where f'(k) is the Fourier transform of fix). The theorem suggests a way to calculate 
spatial derivative fix), as shown in Fig. 3. Firstly, transform f(x) in real space into f(k) 
in reciprocal space; then, multiply f(k) with ik; finally, take the inverse Fourier transform 
(IFT) of f'(k), the spatial derivative f'(x) can be obtained. 

The approach mentioned above has excellent accuracy properties, typically well beyond 
that of standard discretization schemes. In principle, it gives the exact derivative with 



m 


FFT 


f(k) 


xik 


> 


> 



ikxf(k) = f'(k) 



IFFT 


/'(*) 


> 



Figure 3: A possible flow chart for differential operator using the FFT scheme and its inverse (IFFT). 



infinite order accuracy if the function is infinitely differentiable |53|, [54|, [55|, [56|, which is 
another advantage of FFT scheme compared to the FD scheme. In our manuscript, using 
this virtue, the FFT scheme is designed to approximate the true spatial derivatives, as a 
result, to eliminate spurious velocities and to guarantee energy conservation. 

However, the trouble in proceeding in this manner is that, in many cases, it is difficult to 
ensure that infinite differentiability condition is satisfied. For example, the Sod shock tube 
problem j^] contains the shock wave, the rarefaction wave and the contact discontinuity. 
Then the derivative of hydrodynamic quantity, p'(x) or T'(x) has a discontinuity as the 
same character as the square wave (see Fig. 6 for more details). Then the discontinuity will 
induce oscillations, known as the Gibbs phenomenon. The Gibbs phenomenon influences 
the accuracy of the FFT not only in the neighborhood of the point of singularity, but also 
over the entire computational domain. More importantly, sometimes, it will cause numerical 
instability. For example, for the problems shown in Figs. 9 and 17 (in Sec. IV), the 
above approach is unstable due to the Gibbs phenomenon. Recently, there is a trend to use 
smoothing procedures which attenuate higher-order Fourier coefficients to avoid or at least 
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59] . A straightforward and 



to reduce these oscillations (i.e., WFFT method) 
convenient way to attenuate the higher-order Fourier coefficients is to multiply each Fourier 
coefficients by a smoothing factor (filter) <Tfc, such as the Lanzos filter, raised cosine filter, 
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59] 



sharpened raised cosine filter and exponential cutoff filter, as listed in Refs. 

In the present study, based on Taylor series expansion of wave number k, we present a 
way to construct smoothing factors. Firstly, we expand k in Taylor series 
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Figure 4: Smoothing factors <Jk for hi, k2, &3, and k& with N = 128. 

where T(n) = / °° t n ~ 1 e~ t dt is the Gamma function, 0(n) = Mod[— 1 + ra, 2] is the Mod func- 
tion, and e(— 1 + n) is the unit step function. Thus, in order to damp the Gibbs oscillations, 
or in order to filter out more high frequency waves, k may take the form of an appropriately 
truncated Taylor series expansion of sin(/cAx/2). For example, k may take the following 
forms 

sin(fcAx/2) 



ki- 



k 2 =h + 



k 3 =k 2 + 



Ax/2 ' 

sin 3 Ax/2) /6 
Ax/2 ' 

3sin 5 OAx/2)/40 
Ax~/2 ' 



and 



k 4 =k 3 



5sin 7 (fcAx/2)/112 



A 



xj2 



(44) 
(45) 
(46) 

(47) 



where k\ is consistent with the one used in Ref. [60|. Some simple derivations indicate that 
the above approach with k±, k 2 , k%, and k^ has a second-order, fourth-order, sixth-order, and 
eighth-order accuracy in space, respectively (see Appendix for more details). 
Therefore, smoothing factor for k\ takes the following form 

k\ sin(n7r/iV) 



n = -N/2,..., N/2. 



(48) 



k titt/N 

Smoothing factors cifc for k 2 , k%, and k± can be calculated in a similar way and are represented 
in Fig. 4. It is clear that the lower-order smoothing factors k\ and k 2 , filter out more high 
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Figure 5: (Color online) Absolute errors for the computed first-order derivative of the test function f(x) 
sin(x) with the WFFT algorithm. 




Figure 6: (Color online) Comparisons between LB results and exact ones for the one-dimensional modified 
Sod problem at t = 0.2. The simulation results (circles) are from the WT model with the WFFT scheme, 
where the lower-order filter k\ is used. 




X X 

Figure 7: (Color online) Temperature profiles for the modified Sod shock problem obtained from the WFFT 
schemes with k, k\, kz, &3, and k^ in (a). The portion in square of (a) is enlarged in (b) for a closing view. 

frequency waves, and may result in excessively smeared approximations, which are unfaithful 
representations of the truth physics. On the other hand, the higher-order smoothing factors 
&3 and hi, reserve more higher frequency waves, but may not damp the Gibbs phenomenon 
(see Fig. 7 for more details), then cause numerical instability. The smoothing factors should 
survive the dilemma of stability versus accuracy. In other words, they should be minimal 
but make the evolution stable at the same time. 

As a simple test, using the WFFT algorithm, the derivative of a infinite differentiable 
function f(x) = sin(x), is calculated with ki, fc 2 , k 3 , /c 4 , and plotted in Fig. 5. It is clearly 
seen that when k^ is used, the errors reduce to round-off. As another test, the validity of the 
WFFT scheme is verified by the modified Sod shock tube with higher pressure ratio. For 
the problem considered, the initial condition is described by 

(p, u, v, T)\ L = (1.0, 0.0, 0.0, 1.0), x < 0; 

(49) 

(p, u, v, T)\ R = (0.125, 0.0, 0.0, 0.5), x > 0. 

Subscripts "L" and "i?" indicate macroscopic variables at the left and right sides of the 
discontinuity. The size of grid is Ax = Ay = 2.5 x 10~ 3 , time step is At = 10~ 5 , and 
relaxation time is r = 3 x 10~ 4 . Figure 6 shows the computed density, pressure, velocity, and 
temperature profiles at t = 0.2, where the circles are for simulation results and solid lines 
are for analytical solutions. The two sets of results have a satisfying agreement. Figure 7 
shows the temperature profiles obtained from the WFFT schemes with k, k\, &2, k 3 , and k^ 
in (a) and local details of the part near the shock wave in (b). One can see that higher-order 
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Figure 8: (Color online) Variations of total energy Aer(f) for the phase separating process described in Fig. 
2, obtained from the WFFT schemes with k\, &2, ^3, and ki. The ones obtained from the WFFT schemes 
with fc2, &3, and ki are enlarged in the inset. 

filters, such as and ki, have higher accuracy in smooth regions, but cannot refrain the 
Gibbs phenomenon in unsmooth regions effectively. The lower-order filters, although are 
too dissipative, can damp the spurious oscillations to neglectable scale, which are capable of 
shock capturing. Therefore, it should be noted that, for flows without shock waves and/or 
discontinuities, the WFFT scheme with higher-order filter is stable, valid and appropriate. 
While for the compressible flows with shock waves and/or discontinuities, the WFFT scheme 
with lower-order filter is a more appropriate choice. In the present study, we focus on the 
liquid-vapor systems without shock waves and strong discontinuities. Therefore, the WFFT 
schemes with higher-order filters are used. 

For comparisons, we verify the proposed FFT algorithm with the same problem described 
in Fig. 2 and display variations of total energy Aey(t) obtained from WFFT schemes with 
ki, fe, &3, and k± in Fig. 8. It is found that, for each case, Aex{t) oscillates at the beginning, 
then approaches a nearly constant value. Behaviors of A&r(t) can be interpreted as follows. 
At the beginning of phase separation, the fluids separate spontaneously into small regions 
with higher and lower densities, and more liquid- vapor interfaces appear. As a result, spacial 
discretization errors in Eq. (10) induced by interfaces arrive at their maxima that account for 
the initial oscillations. As time evolves, under the action of surface tension, the total liquid- 
vapor interface length decreases due to the mergence of small domains, then the discretization 




errors decrease. With the increase of precision, variations of total energy Aetoi{t) decreases. 
We can, therefore, come to the conclusion that WFFT scheme with higher-order filter has 
more advantage in guaranteeing energy conservation than the one with lower-order filter and 
other FD schemes used above. 

4. Simulation results and analysis 

In this section, two kinds of typical benchmarks are performed to validate the physical 
properties of the thermal multiphase model and the newly proposed algorithm. The first one 
is related to a planar interface. The second one is related to a circular interface. 

4-1- Coexistence curve, spurious velocities and interfacial width 

To check if the thermal LB multiphase model can correctly reproduce the equilibrium 
thermodynamics of the system and the numerical accuracy of the new scheme, a series of 
simulations about the liquid- vapor interface at different temperatures were performed. Unless 
otherwise stated, the WFFT scheme with is used throughout our simulations. 

Simulations were carried out over a 512 x 4 domain with PBC in both directions. The 
initial conditions are set as 

' (p, u x , u y , T)\ L = ( Pv , 0.0, 0.0, 0.97), x<N x /4; 
< (j>,u x ,u v ,T)\ M = (p l , 0.0,0.0, 0.97), N x /A < x < 3iV x ./4; (50) 
k 0, u x , u y , T)\ R = ( Pv , 0.0, 0.0, 0.97), 3N X /A < x, 

where p v = 0.80 and pi = 1.20 are the theoretical values at T = 0.99. Parameters are set to 
be t = 10~ 2 , K = and others are unchanged. The initial temperature is set to be 0.97, but 
dropping by 0.02, when the equilibrium state of the system is achieved. Simulations were 
then run until the temperature had reduced to 0.87. 

In Fig. 9, the liquid- vapor coexistence curves from LB simulations using various numerical 
schemes at different temperatures are compared to the theoretical predictions from Maxwell 
construction. One can see that when using the WFFT and the NPS schemes, the results 
are closer to the theoretical phase diagram. Nevertheless, when the temperature is lower 
than 0.87, the NPS scheme becomes unstable. Physically, this is owing to the sharp interface 
when K = (see Fig. 11) and the nature of the VDW EOS, since large density ratio 





Figure 10: (Color online) Velocity profiles obtained using various schemes at T = 0.93 with K = 0. 
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Figure 11: (Color online) Density profiles obtained using different schemes at T = 0.93 with K = 0. 
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Figure 12: (Color online) Velocity profiles obtained using various schemes at T = 0.93 with K = 5 X 10~ 6 . 

occurs when the temperature is much lower than the critical one. From another perspective, 
it demonstrates that the WFFT algorithm has a better numerical stability for this test. 
Results from the MLW2 and the FL schemes deviate remarkably from the theoretical values, 
especially for the vapor branch at lower temperatures. Besides physical reasons listed above, 
numerical accuracy of these two schemes is also an important factor. 

The velocity and density profiles at T = 0.93 are shown in Fig. 10 and Fig. 11, respec- 
tively. As one can see in Fig. 10, for all schemes, spurious velocities exist and reach their 
maxima near the interface regions. However, the maximum spurious velocities obtained from 
different schemes are greatly different. For the MLW2 and the FL schemes, the maxima of u s 
are on the order of 10 -4 and 10 -2 , respectively. A significant reduction of the u s is achieved 
by using the NPS scheme, decreasing the maximum to about 5 x 10~ 7 . Through the usage 
of the WFFT algorithm, u s is further reduced by an order of magnitude compared to NPS 
scheme. Density profiles in Fig. 11(a) indicate that spurious interfaces (scatter symbols near 
the interfaces) have been produced when using the MLW2 scheme or the FL scheme because 
of excess numerical diffusion, which does not provide us a clear picture of phase separation, 
especially when the temperature is close to the critical value. This feature is not present 
when using the NPS scheme or the WFFT scheme (see Fig. lib). 

For all numerical schemes, the strength of surface tension plays an important role in 
reducing spurious velocities, as shown in Fig. 12. In the case of K = 5 x 10~ 6 , the amplitudes 
of u s can be reduced by a factor of approximately 10 with respect to the case of K = 0. 
Subsequent simulations indicate that w^ax wm decrease to 10~ 12 when K increases to 10~ 5 . 
This is due to the existence of a wider interface and a smaller density gradient in the interface 
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Figure 13: (Color online) Comparisons of coexistence curves from LB simulations and Maxwell construction. 
Here K is set to be 5 x 1CT 6 . 
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Figure 14: (Color online) Velocity profiles obtained using various schemes at T = 0.85 with K = 5 X 10 



region when K increases. With the decrease of spurious velocities, a more accurate phase 
diagram, especially in the vapor branch, is also achieved (see Fig. 13), even in the MLW2 
case and the FL case. 

Besides the strength of surface tension, temperature is another key factor affecting spuri- 
ous velocities, as displayed in Fig. 14. For all numerical schemes, a lower temperature makes 
larger spurious velocities. It is worth mentioning that, at the same temperature, the WFFT 
algorithm also allows to reduce u s of about an order of magnitude compared to NPS scheme. 

In our simulations, so far, we have not discussed in detail the width of interface. According 
to the VDW theory, the interface width I, can be determined by numerically solving the 
following integral for a planar interface 6l|, [62] 

1 f p *^ dp* 

~V2(a/Ky/ 2 J p * (x * o) [$*(p*)-$*(p*)]i/2' 



I = X — Xq 



(51) 




Figure 15: (Color online) Equilibrium density profiles across liquid-vapor interface at T = 0.93 with K = 
5 x 10~ 6 . 

where p* = pb, T* = bT/a and, 

$* = p *£ - p*T*[\n(l/p* - 1) + 1] - p*\ (52) 

t = T* \n(l/p* s - 1) - p* s T*/(l - pi) + 2 P :,s = v,l (53) 

with a — 9/8, 6 = 1/3 in this model. Note that the solution of the above equations gives the 
exact density profile for a planar interface for any value of T. Equilibrium density profiles 
across the liquid-vapor interface at T = 0.93 from LB simulations versus results from VDW 
theory are shown in Fig. 15. It is clear that although the liquid and the vapor densities 
calculated from the MLW2 and the FL schemes coincide with the theoretical ones, neither the 
MLW2 nor the FL scheme produces the correct interface profile. The wider interface in these 
two cases is due to the excess numerical diffusion, as shown in Fig. 11a. The NPS approach 
leads to a small deviation from the VDW theory, while the WFFT scheme presents a perfect 
consistency with the theoretical solution. In Fig. 16 we display density profiles obtained 
from WFFT algorithm with K\ — 5 x 10" 6 and K2 = 1CT 5 , respectively. As expected, the 
interface becomes wider as K increases. A wider interface decreases the density gradient in 
the interface region and helps to stabilize the liquid-vapor system at lower temperature. 

4-2. Laplace's law 

In this subsection, we will look at the dynamics of the relaxation of a deformed droplet 
driven by surface tension, and investigate the magnitude, as well as the spatial extent of 




Figure 17: (Color online) Evolution of a droplet from triangle to circle, where t = in (a), t = 0.1 in (b), 
t = 1.5 in (c), and t = 5.0 in (d). 




Figure 18: (Color online) Velocity fields at t = 30.0 from the NPS scheme in (a) and the WFFT scheme with 
&4 in (b). 

the spurious currents for the circular interface. Initially, an equilateral triangular droplet 
with an initial sharp interface is placed at the center of the computational domain with 
N x = N y = 128 lattice units. Initial conditions are given by 

(p, u x , u y , T)\ in = (1.46, 0.0, 0.0, 0.95), 

(54) 

(p, u x , u y , T)U = (0.58, 0.0, 0.0, 0.95), 

where subscripts 'in' and 'out' indicate macroscopic variables inside and outside the liquid 
drop, respectively. PBC are employed on both the vertical and horizontal boundaries. The 
surface tension parameter is K = 10 -5 , leaving the others unchanged. After 3 x 10 6 time 
steps, the system reaches equilibrium. Contour plots of the fluid density at four representative 
times are shown in Fig. 17. It is clearly seen that due to the effects of surface tension, the 
droplet relaxes to a circle slowly. 

The velocity fields at t = 30.0 obtained from the NPS scheme and the WFFT scheme 
with &4 are plotted in Fig. 18. To illustrate the structure of the velocity field clearly, the 
lengths of the velocity vectors are multiplied by 5 x 10 5 . To be seen is that spurious currents 
exist in each case and are roughly aligned in the direction normally to the interface and 
rapidly disappear away from the interface. However, the magnitude of the spurious currents 
are significantly reduced as the WFFT approach is used. Figure 19 shows temporal evolution 
of the maximum velocity u s max with the second-order, the fourth-order, the sixth-order, and 




Figure 19: (Color online) Evolution of the maximum velocity u s max versus time t with the NPS scheme and 
the WFFT schemes with k%, k%, fe, and Ha- 

the eighth-order FFT schemes, and the NPS scheme. We can see that, in each case, u s max 
decreases, and tends to nearly a constant when t > 20. More importantly, with the increase 
of precision, u s max decreases. There is a decrease of a factor 22 for the velocities when using 
the WFFT scheme with respect to the NPS scheme. 

The density and the pressure profiles along the center line of the droplet are plotted in 
Fig. 20. In the inner and outer of the droplet, pressures P; n and P ou t are two constants 
and a rapid change occurs across the interface. The difference between the two constants is 
usually used to compute the surface tension for a given K. For this purpose, we introduce 
the Laplace law which states 



a = rAP = r{P in - P ( 



out ) 



(55) 



where Pi n is the mean pressure inside the droplet averaged over all points of Ti n < r/2 from 
the droplet center and P out is the external mean pressure averaged over all the points of 



r out > 3r/2. In this way, only the particles 
tension can also be computed in such a way 

a = K 



ar from the interfaces are considered. Surface 
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(56) 



In order to test these relations, a series of simulations with sides ranging from 48 to 81 
are run with three different surface tension parameters K\ = 10~ 5 , A'2 = 7.5 x 10 -6 and 
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Figure 20: (Color online) Density and pressure profiles along the center line of a droplet. 
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Figure 21: (Color online) Laplace law tests for three different surface tension parameters. The scattered 
symbols are for simulation results and the dashed lines are linear fits of the scattered symbols. 

K% = 5 x 1CT 6 . In Fig. 21, we present a plot of AP versus 1/r at T = 0.95 and linear relation 
is well satisfied. By measuring the slope, the surface tension is found to be o\ = 1.98 x 1CT 4 , 
cr 2 = 1.77 x 10 -4 , and 03 = 1.34 x 10 -4 , which are in excellent agreement with the theoretical 
values, obtained from Eq. (55), o x = 1.96 x 1(T 4 , a 2 = 1.72 x 10" 4 , and a 3 = 1.40 x 10" 4 . We 
mention that the relative error of surface tension e = \<Jlb — <5exact \ Inexact increases with the 
decrease of surface tension parameter. There are two reasons accounting for this behavior. 
Firstly, a larger K will cause a larger surface tension and a larger pressure difference. This 
helps to measure the pressure difference with high accuracy. Secondly, when K decreases 
larger spurious velocities will be produced and larger pressure oscillations will be induced. 





5. Conclusions and discussions 

In this paper, a thermal LB model for liquid-vapor system is developed. The present 
model experienced mainly three stages. It was originally composed by WT for ideal gas, 
then developed by GLS by adding an interparticle force term. Here we propose to use the 
WFFT scheme to calculate both the convection term and the external force term. The usage 
of the WFFT scheme is detailed and analyzed. It is found that the lower-order filters, with 
better numerical stability and lower accuracy, can effectively reduce the Gibbs phenomenon 
at discontinuity, while the higher-order filters help the scheme to maintain high resolution in 
smooth regions. One can choose appropriate filter according to the specific problem. With 
the higher-order WFFT algorithm, one can better control the non-conservation problem of 
total energy due to spatiotemporal disterizations. The model has been successfully applied 
to the calculation of interfacial properties of liquid-vapor systems. Very sharp interfaces 
can be achieved. By adopting the new model the magnitude of spurious currents can be 
greatly reduced. As a result, the phase diagram of the liquid- vapor system obtained from 
simulations are more consistent with that from theoretical calculation. The accuracy of the 
simulation results is also verified by the Laplace law. 

Besides the numerical effects, both the surface tension and temperature have also signifi- 
cant influences on the spurious velocities. A stronger surface tension and/or a higher temper- 
ature can decrease the density gradient near the interfaces and stabilize the simulations. The 
analysis presented in this work provides an convenient way of extending the WFFT approach 
to multiphase LB models and to numerical solving partial differential equations. In further 
studies we will increase the depth of separation, which the model can undergo and investigate 
the similarities and differences between thermal and isothermal phase separations. 
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Appendix A. Appendix: numerical accuracies of operators fci, k%, ks, and k± 

Substituting Eq. (S3]) into Eq. flU, the RHS of Eq. can be expressed as 

ifcx/>) = ^[sin(fcAx/2) x f(k) 

+ 1 sin 3 (k Ax/2) x f(k) 
6 

+ ^ sin 5 (k Ax/2) x f(k) 

+ ±-sm 7 (kAx/2)xf(k)---}, (A.l) 
Taking IFT of the RHS of the first line of Eq. ([Oil gives 

IFT[i/ci x J{k)\ = - ^ e i ^ x^-sin(A;Ax/2)x/(A ; ) 

n=-N/2 ' 
1 N/ 2 - 1 ifcAx/2 _ -ifcAx/2 _ 

L ^ Ax w 

n=-N/2 

f(.r, ■ Ar-2) /(,■, -A.r 2) 
Ax 

= f'ixj) ■ ^A.Sf'Uj) ■ ... (A.2) 

It is clear that the FFT scheme with operator ki has a second-order accuracy in space. 
In a similar way, we have 

IFT[ifc 2 x /(*)] = f\x 3 ) + J^A* 4 /^,) + -> (A.3) 

IFT[ifc 3 x f{k)\ = f{ Xj ) + ^l^Ax 6 / (7) (^) + (A.4) 
IFT[ifc 4 x f(k)] = f( Xj ) + 9289 1 7280 Aa;8 / (9) (^) + (A.5) 



where f^(xj), f( 7 \xj), and f^ 9 \xj) represent the fifth order, the seventh order, and the 
ninth order derivatives, respectively. Therefore, the WFFT approach with ki, k 2 , k%, and 
k& has a second-order, fourth-order, sixth-order, and eighth-order accuracy in space, respec- 
tively. 

From another perspective, it should be noted that, the FFT scheme is not a local scheme 



or, in other words, k is not a local operator 55 



since each FFT coefficient is determined 



by all the grid point values of f(xj), as shown in Eqs. (l401l4Tl) . Therefore, the FFT scheme is 
not a finite-point formula, like the second-order FD is a 3-point formula, or the fourth order 
expression, is a 5-point formula; rather, the FFT scheme is iV-point formulas. But there 
are important reasons for expressing derivatives as local operators. In a continuous space, 
the derivative of a function is defined locally. Hence, when modeling a continuous system 
with a discrete system, it is desirable to retain the local character of the derivative. This 



can be especially true near boundaries or marked internal inhomogeneities 54|. From the 
above derivations, we find that the FFT scheme with ki corresponds a 3-point FD scheme. 
Therefore, from the point of numerical analysis, kx, k 2 , k^, and k^ can maintain the local 
characteristic of k in some extent. Hence, errors arising from the discontinuity are also 
localized and the accuracy away from the discontinuity can be ensured. 
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